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The semiclassical theory for billiards with mixed boundary conditions is devel- 
oped and explicit expressions for the smooth and the oscillatory parts of the spectral 
density are derived. The parametric dependence of the spectrum on the boundary 
condition is shown to be a very useful diagnostic tool in the semiclassical analysis 
of the spectrum of billiards. It is also used to check in detail some recently pro- 
posed parametric spectral statistics. The methods are illustrated in the analysis of 
the spectrum of the Sinai billiard and its parametric dependence on the boundary 
condition on the dispersing arc. 

1 Introduction and Statement of the Main Results 

A classical billiard is completely defined once its boundary E is prescribed. Then, the 
requirement that the particles (rays) reflect specularly determines the dynamics. To address 
the quantum (wave) analogue, it is necessary to solve the Schrodinger equation, which 
reduces for billiards to the Helmholtz equation 



k is the wave number under consideration and we use natural units where h = 2m = 1 and 



In wave dynamics, one needs to supply an additional piece of information, that is, a 
condition which the wave function has to satisfy on the boundary. In quantum mechanical 
applications, one often considers the boundary as high potential wall whose height tends 
to infinity and then the boundary condition is the Dirichlet condition 
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Another common boundary condition is the Neumann boundary condition 
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where dn stands for the normal derivative, with the normal pointing outside. This boundary 
condition appears most naturally for the pressure field in acoustics. One can generalize the 
above boundary conditions by requiring 

K{f)il^{f) + dni^{f) = {) , f ^Y. . (4) 

This form leaves the problem self-adjoint when n is real. The Dirichlet boundary condition 
is recovered by setting n{r )~^ = 0. We shall restrict the discussion to non-negative n{r ) 
functions. If n{r ) is negative, the Helmholtz equation can have solutions for imaginary 
values of the wave number k. Such solutions are unphysical, since they correspond to 
steady states despite of a dispersion of energy in the medium that takes place when the 
imaginary part of k is different from zero [I]. Furthermore, we shall not deal with the 
most general positive n{r) but restrict our attention to piecewise constant functions on the 
boundary. In some applications it is convenient to express eq. (4) in terms of a positive 
parameter h and a (piecewise constant) mixing angle a so that the boundary condition 
reads 

h cos ail){r) -\- sin a dnt/j[f) = , r G S . (5) 

This form interpolates conveniently between the Dirichlet and the Neumann boundary 
conditions, in terms of a which is limited to the interval < a < I-. 

At this point we would like to make the following observation. In the semiclassical 
domain, the dominant term in the boundary condition (4) is dnt/j[f) which is of order k. 
Hence, for a fixed k, and in the semiclassical limit [k — )• oo), the spectrum will always tend 
to the Neumann case. A more proper definition of the semiclassical limit is one, that allows 
also in the limit A; — )• oo an interpolation between the Dirichlet and Neumann boundary 
conditions. This can be achieved by considering k and n/k as independent parameters 
when performing this limit. To emphasize this point, we shall always quote the results as 
functions of these two parameters. 

Mixed boundary conditions are not encountered in physical applications as often as 
Dirichlet and Neumann boundary conditions. A perturbative treatment of the Helmholtz 
equation with respect to changes in boundary conditions can be found in [2], but the semi- 
classical quantization of billiards with mixed boundary conditions was rarely discussed 
previously. The only exception we found was the work of Balian and Block, who encoun- 
tered these boundary conditions in a nuclear physics context [3]. 

Our interest in the mixed boundary conditions stems from the observation that one can 
use the additional freedom provided by the function K,[f) as a powerful tool in the analysis 
of various aspects of the semiclassical quantization of billiards. To explain this rather non- 
conventional approach, we shall have to start by quoting the results of the semiclassical 
derivation for two-dimensional billiard systems which is presented in the next two sections. 

We consider the spectral density in terms of the wave number k. For positive k and 
knin) = +^En{n) it is given by 

OO 

d{k; k) = ^8{k - kn{n)) = d{k; k) + dosc{k; k) , (6) 

which depends parametrically on the boundary conditions. The semiclassical theory uses 
separate methods to evaluate the smooth and the oscillatory components of the spectral 
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density. For the smooth spectral density we get (see sections 2.1, 2.2, 3.1 and appendix A 
for the derivation) 
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This is the generahzed Weyl formula. As expected, the leading term involving the area 
of the billiard A is independent on the boundary condition. The higher order corrections, 
starting from the term containing the circumference Z, depend on /€ in a way which in- 
terpolates between the known expressions for the Dirichlet and the Neumann boundary 
conditions. The third term is the curvature contribution which contains an integral over the 
curvature l/i?(s), and dc{k; k.) denotes contributions from corners. For a 90° corner with 
mixed boundary conditions on one side and Dirichlet or Neumann boundary conditions on 
the other side it is given by 



where the negative sign corresponds to the Dirichlet case. 

The semiclassical expression for dosc{k] n) is written as a sum of oscillatory terms whose 
periods correspond to the lengths of classical periodic manifolds. The semiclassical treat- 
ment distinguishes between contributions of unstable, isolated periodic orbits and contri- 
butions of continuous manifolds of neutral periodic orbits. In the standard theory, the 
former are given by the Gutzwiller trace formula [4, 5] and the latter were first derived by 
Berry and Tabor [6, 7]. To leading order, the introduction of mixed boundary conditions 
does not affect the amplitude of the oscillating terms. The phase of each term is changed 
relative to the case with Dirichlet boundary conditions, by the addition of the n dependent 
phase 



where the periodic orbit bounces n times off the boundary at the points {ri}i=i^...^„ and 
kcosO^ is the component of the momentum normal to the boundary at the i'th bouncing 
point. This expression is the same whether the periodic orbit is an isolated unstable 
periodic orbit (see equations (77) and (81) ) or belongs to a manifold of neutral periodic 
orbits (see equations (37) and (47) ). When K,[r) = 0, one recovers the well known result 
that there is a phase difference of mr between the contributions of periodic manifolds in the 
Dirichlet and Neumann cases. An important feature of (9) is that except for the limiting 
Dirichlet and Neumann cases the phase depends on the normal momenta at the bouncing 
points. We would like to emphasize that the underlying classical dynamics is indifferent 
to the boundary conditions. In other words, by changing the boundary condition, one can 
alter the quantum spectrum, without affecting the periodic orbits. 
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We shall explain now how the parametric dependence discussed above can be used as a 
diagnostic tool in the semiclassical analysis of spectra. Suppose we would like to isolate the 
contribution of periodic orbits which are confined to a certain segment of the boundary. We 
may calculate the spectrum for two K,[r) functions which are different only in the desired 
segment. The difference in the oscillatory parts of the corresponding spectral densities will 
then depend only on trajectories which bounce off the domain of interest at least once. This 
is of particular importance when the boundary under consideration permits the existence of 
continuous families of marginally stable periodic orbits such as the "bouncing ball" orbits 
in the Sinai and the stadium billiards [8, 9]. Their effect is best seen when the spectral 
density is Fourier transformed to give a length spectrum. Here, each family contributes a 
peak at the lengths of the bouncing ball orbit and its repetitions. If the original spectrum 
is considered in the interval [0, k^ax]} the intensity of the peaks at the bouncing ball lengths 
for a J-dimensional billiard scales like k^~^'^/''^ relative to the peaks appearing at the lengths 
of the unstable periodic orbits. Moreover, as d increases the number of possible bouncing 
ball families proliferates and they can fill subspaces of dimensions ranging between two and 
d in configuration space. Thus, the bouncing ball contributions may dominate the length 
spectrum to the extent that contributions from the unstable periodic orbits can hardly 
be resolved. The variation of the boundary conditions along the sections of the boundary 
which are visited by the unstable periodic orbits exclusively, enables us to isolate their 
contribution from the non-generic features which are due to the bouncing ball families. We 
shall show the power of this method in section 4 where we discuss in detail the spectrum 
of the Sinai billiard in the plane. 

Another interesting application is the possibility to study the parametric spectral statis- 
tics, along the lines which were introduced in [10, II, 12, 13] for parameter dependent 
Hamiltonian systems. The system that we study is different from the standard one, since 
the underlying classical dynamics is independent on the parameter a, so that the classical 
periodic orbits are the same for all parameter values. The only change in the semiclassical 
expression is due to the phase factor (9). In particular, we shall study in section 5 the dis- 
tributions of spectral "velocities" dEn{a) / da and "curvatures" d^En{a)/da^ for the Sinai 
billiard. Parametric statistics of this kind were previously performed for parameter depen- 
dent Hamiltonian systems and universal distributions for chaotic systems were suggested, 
based on random matrix models. 

Much of the present work was dedicated to the derivation of the semiclassical spectral 
density for the mixed boundary condition. This requires some special care, as was first 
noted by Balian and Bloch in their work on the smooth spectral density for billiards in 
three dimensions with mixed boundary conditions [3]. In their analysis, Balian and Bloch 
explain very clearly the main difficulty which is encountered when one addresses general 
boundary conditions rather than the Dirichlet or Neumann conditions. Their starting 
point is an integral equation for the Green function, with a kernel which involves the free 
Green function and its normal derivative on the boundary. A Born expansion then leads 
to a multiple refiection series for the Green function. Balian and Bloch showed that the 
kernel of the integral operator is regular for Dirichlet and Neumann boundary conditions. 
However, for any intermediate case the kernel is singular (but still integrable). Because of 
this singular behavior, any perturbative expansion or semiclassical approximation, which 
can be used for the two extreme boundary conditions cannot be justified for any interme- 
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diate case. Balian and Bloch observe that "It is therefore not possible to use directly the 
integral equation for deriving a perturbation expansion. . . For instance, even for a plane, 
the expansion should be resummed to all orders in k." Balian and Bloch did not provide 
the corresponding treatment for the oscillatory component of the spectral density, nor did 
they compute the smooth spectral density for two-dimensional billiards. Moreover, their 
method excludes the treatment of corners. The following section will be devoted to a dis- 
cussion of these points. We shall use a variety of methods to overcome the difficulty which 
was identified by Balian and Bloch. Sections 2 and 3 are rather technical. The reader 
who is not interested in the details of the derivation can skip them and, equipped with 
equations (7,9) he/she can go directly to the applications. 

The paper is organized in the following way. In the next two sections we shall present the 
semiclassical theory for the spectral density. We start with integrable billiards, the circle 
and the rectangle, and then proceed to more general shapes which correspond to chaotic 
classical billiards. The numerical demonstrations and checks will be carried out for the 
Sinai billiard and therefore we shall give an independent derivation of the spectral density 
using a method which is based on the KKR technique as used by Berry [8]. In all cases, we 
shall discuss separately the smooth and the oscillatory parts of the spectral density. For 
the former, we shall use a variety of methods, extending the Stewartson and Waechter [14] 
calculations for the circle and the Balian and Bloch approach for smooth billiards. The 
resulting expression (7) for the smooth density includes the modifications due to the mixed 
boundary conditions for the length and curvature terms as well as a 90° corner term. To 
calculate the modification of the oscillating part of the spectral density we shall mostly 
use the scattering approach [15, 16, 17]. Section 4 will be devoted to the application of 
the mixed boundary conditions as a diagnostic tool. We shall develop a method for the 
elimination of all the structures which are due to the bouncing ball families and show how 
it works in practice for the Sinai billiard in two dimensions. In section 5 we shall discuss 
spectral statistics that depend on the sensitivity of the energy levels to a change of an 
external parameter. In particular, we compare numerical results to distributions that were 
obtained from random matrix theory, and we discuss deviations from these distributions 
that are due to the existence of families of bouncing ball orbits. Again, the system to 
be analyzed is the Sinai billiard and the boundary condition on the dispersing arc will 
be varied. The last section contains some discussion about the physical interpretation of 
the mixed boundary conditions and the behavior of the boundary phase for long unstable 
periodic orbits. 

2 The Semiclassical Spectral Density — Integrable 
Systems 

The semiclassical theory relates the oscillatory part of the spectral density to the peri- 
odic manifolds of the underlying classical dynamics. This is why it is necessary to treat 
separately billiards whose classical analogues are integrable or chaotic. We shall follow 
this route in this paper and present now the theory for two integrable billiards - the circle 
and the rectangle, deferring the treatment of chaotic billiards to the next section. The 
semiclassical expressions for the smooth components of the spectral density do not depend 
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on the detailed classical dynamics. Because of the relative simplicity of the circle and the 
rectangle billiards one can use special methods to derive the smooth spectral densities in 
these cases. We quote them here since they provide both physical insight and independent 
checks for the general theory presented in the next section. Moreover, the treatment of 
corners will not be done by the method of Balian and Bloch and we rely on the analysis of 
the integrable problems for this information. 

2.1 The Circle Billiard 

We start with the circle billiard with mixed boundary conditions which is integrable if K,[f) 
is constant. In this case the semiclassical approximation can be derived from the exact 
solutions for the Green function and the scattering matrix. Although the circle billiard 
is a special system, the results already show how the semiclassical trace formula will be 
modified in the general case. 

The exact solutions of the Helmholtz equation for the circle billiard with a constant 
K,[f) are given in polar coordinates by 

i^iAr, 0) = ci,m Mhmr) e''' , / G Z , m G Z\{0} , (10) 

where are normalization constants and Ji{z) are Bessel functions of the first kind. The 
wave numbers ki^m are determined by the boundary conditions (4). They have to satisfy 

KJl{kl,mR) + h,m J'i{h,mR) = . (11) 

For < /€ < oo all solutions of eq. (11) are real [18] and the positive solutions ki^ra lie in 
between those for Neumann and Dirichlet boundary conditions: kf^^ < ki^m < kf^ , m > 
0. There exists an equal number of negative solutions ki^^m = —ki^m which correspond 
to the same quantum state. Let us discuss briefiy, what happens if k, becomes negative. 
For every /, the first positive solution ki^i of eq. (11) decreases as k, is decreased and it 
coincides with zero and with the first negative solution ki^^i if k, = —1. If /€ is decreased 
further, then these two solutions become a pair of complex conjugate imaginary numbers. 
For / = these solutions are imaginary for any negative value of k. Imaginary solutions of 
this kind exist also for other billiards and for that reason our considerations are restricted 
to non-negative values of k, [1]. 

The smooth part of the spectral density 

In order to obtain the smooth parts of the spectral density and of the spectral staircase 
for the circle billiard we apply the method of Stewartson and Waechter [14] who derived 
the corresponding results for Dirichlet boundary conditions. The method is based on the 
fact that the exact Green function of the circular billiard can be written in a closed form. 

The starting point is the Green function for the heat diffusion equation which satisfies 
the inhomogeneous differential equation 

{V^ -s^)G{f,f\s^) = 6{f-f') (12) 

with mixed boundary conditions on a circle of radius R: 

{K + dn)G{f,f',s^) = , \f\ = R. (13) 
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This Green function is directly related to the Green function of the Helmholtz equation by 
G{r,r\s^) = G{r,r\k'^)\k=is and to the spectral density by [k > 0) 



2k 

dik) = lim Tr Im G{f,f',s^ — ii/M^- 

= lim / d^f [Im Gfr, r", — if ) 



-ik 



|s= — «/vJr '= r •} 



(14) 



where the integration extends over the domain of the circle. In the following an asymptotic 
expansion of the Green function for large s will be derived and it yields an asymptotic 
expansion of the smooth part of the spectral density for large k. A detailed discussion of 
this point is given in [19]. 

The solution to equations (12) and (13) is obtained by taking the free Green function 



Go{f,f\s^) = -—Ko{s\f - r 
Ztt 



1 oo 

— J2 In{sr<)K.,{sry)cos[n{e-e')], (15) 



which satisfies (12) and adding to it solutions of the corresponding homogeneous differential 
equation in order to satisfy the boundary conditions. In eq. (15) (r, ^) and {r\6') are the 
polar coordinates of r*and r*', r< (r>) is the smaller (greater) of r and r' . In{z) and Kn{z) 
are modified Bessel functions [20]. The Green function G{r,r\s^) can then be written as 



G{f, r s^) 



I OO 

— In{sr<)[K^{sry) + aJ^{sry)]cos[n{e-e')], (16) 



and the coefficients a„ follow from the boundary condition (13) as 

"""" Kh{sR) + sI'^{sR) 



(17) 



In order to obtain the spectral density one has to evaluate the trace of the Green function. 
However, since the real part of this trace is divergent one considers instead 



R 



2tt 



dr r I d9 [G{f,f' , s'^) — Go{f,f' , s'^)]f 
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which is finite. Evaluating the two integrals results in 
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:i8) 



(19) 



where 
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)In{sR)Kn{sR) - l'SsR)K{sR) 



sRInisR) 



K, -\- S 



sH) 



K„(sR) 



, I'JsR) 



(20) 



The function /(n, s) is even in the first argument for integer n. In order to derive the 
asymptotic behavior of K{s'^) for large s, the sum in eq. (19) can be replaced by the 
integral 



K{s') = R' / dz./(z.,3) 
Jo 



(21] 
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since the correction terms are exponentially small for large s [14]. One proceeds now in 
the following way. The definition (20) of /(n, s) consists of two factors. The first one 
is the expression for /(n, s) in the case of Dirichlet boundary conditions on the circle 
and for i? = 1 it is identical to the definition of Stewartson and Waechter [14]. The 
Bessel functions in this term are replaced by their uniform approximation [20] and after 
substituting p/[sR) by the term is expanded in powers of l/[sR). The same is done 
with the quotients K'^[sR) / Kn{sR) and I'^[sR) / In{sR) in the second factor of definition 
(20). Taking into account terms which contribute to the first two leading terms of K{s'^) 
yields 



K( 



R' 



d/i 



1 



2si?(l + /i^ 



SR{1 + /i2)3/2 



/i2 - [2si?(l + 



f + yrT7?- [23i?(l + /i2)]- 



(22) 

In this expression the variable s appears in two different combinations, as sR or s/k. In 
keeping with our remarks in the introduction we expand the integrand in powers of l/[sR) 
keeping the ratio k/s fixed. This leads to a resummation of d[k) to all orders in /€, which 
is identical to the resummation performed by the method of Balian and Bloch [3] which is 
reported in section 3.1. The result is 



K(s') = R' d/i 
Jo 



+ + 1 



2sR (1 + /i2) (^ + yrTT?) 2(3i?)2 (1 + /i2)5/2 + yrT7?)2 

From this the first asymptotic terms of the smoothed level density d[k) follow as 
d{k) -- 



(23) 



kR^ 
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TT 


kR^ 
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~ ~2 



l + (f 



+ 



1 i+nt 



i+(f 



3/2 



2k 



l + (f 



3/2 



(24) 



where the area term kR^ /2 is obtained from the contribution of the free Green function to 
the trace of G. The integral over fi has been performed using MAPLE. Integrating eq. (24) 
over k yields 



N{k) 



{kRf 



kR 




1 

+ 6 



l + (f 



1 



1 + 



(25) 



In principle, the integration of dik) over k could lead to a K-dependent constant in the 
asymptotic expression of N[k). One way to determine the constant term of the asymptotic 
N[k) is to expand —K[s'^) in all orders of 1/s (now taking into account also the s//€-terms). 
Then the coefficient of the l/s^-term is identical to the constant term in the asymptotic 
series for N[k) [21], if expanded in all powers of 1/k. Doing this, one obtains in both 
cases 1/6 — k i?, and thus the constant term in eq. (25) is correct. Eq. (25) is identical 
to the formula for general shapes that will be derived in section 3, if the three terms are 
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expressed by the area A = vri?^, the perimeter L = 2ttR and the integrated curvature 
Jds / R{s) = 27r, respectively. 

In appendix A the corresponding result for a semi-circle billiard is derived with mixed 
boundary conditions on the circle and Dirichlet or Neumann boundary conditions on the 
diameter. In this way the contribution of a 90° corner is obtained that will be used in 
applications in section 3.3. 

In figure 1 we show a comparison of the three contributions to the asymptotic expression 
(25) for the mean spectral staircase N[k) with numerical data for the circle billiard with 
K, = 150. All 22 540 eigenvalues below wave number k = 300 have been included. In the 
upper figure the area term is compared to the numerical result for N[k) and on this scale 
one cannot see any difference. The intermediate figure shows the length term in comparison 
with the difference between the spectral staircase and the area term. The numerical data 
have been smoothed over the range of approximately 100 levels. Both curves are in good 
agreement over the whole range of k. In order to demonstrate the effect of the resummation 
in /€, we also plotted the length term for Neumann boundary conditions (dashed), since this 
is the next-to-leading asymptotic term for all values of /€, if no resummation is performed. 
The dashed and the full line become parallel for large values of A;, but one can see that the 
resummation is necessary in order to have a good description of the spectrum also for lower 
values of k. In the lower figure the curvature term is compared to the spectral staircase from 
which the area and length terms have been subtracted. In this case the numerical curve 
has to be smoothed strongly (over the range of approximately 1000 levels) since it has large 
fiuctuations about its mean (note the different scales in the three figures). Nevertheless, 
one gets good agreement with the theoretical curve. The deviations near k = 300 and 
k = (also in the intermediate figure) are due to edge effects of the smoothing procedure. 

The oscillatory part of the spectral density 

We derive the oscillatory part of the spectral density by applying a scattering formalism 
and making use of the fact that there is a direct relation between the inside and outside 
problem of a billiard system. 

A solution of the scattering problem for a compact billiard system can outside a circle 
that contains the billiard completely be written as 

oo 

,l;{r,e) = Hf{kr)t'e^"+ ^ Si,i,{k)H+{kr)t''e^''' , (26) 

/' = -oo 

where Si^ii{k) are the elements of the on-shell scattering matrix in the angular momentum 
representation, and H^[z) = Ji{z) it iYi{z) are Hankel functions. For the circle (with 
radius i?), the requirement that t/j[r^6) satisfies the boundary conditions (4) with constant 
K, gives the scattering matrix explicitly. Due to the rotational invariance of the system, the 
scattering matrix is diagonal, and its diagonal elements are given by 

KHi{kn) + k Hi [ kK) 

where ^i{k) are the scattering phases. From this equation, and from the eigenvalue equation 
for the inside problem (11), it follows that 

Si{k) = l ^ A; G n = 1,2,...} (28) 
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300 



Figure 1: The area, length and curvature term ol the mean spectral staircase N[k) (full 
lines) for the circular billiard with k, = 150 in comparison with the spectral staircase N[k) 
from which the respectively stronger asymptotic terms have been subtracted (dotted lines). 
In the second and third figure the dotted lines have been obtained by smoothing the original 
curves. The dashed line shows the length term for Neumann boundary conditions. 
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for all / G Z. This is a direct relation between the inside problem and the scattering 
problem, and is a particular example of the inside-outside duality for billiard systems [22], 
which has been proved in a general form recently [17]. 

The semiclassical approximation for the S'-matrix is obtained by using the Debye ap- 
proximation for the Hankel functions (or equivalently by a WKB approximation). This 
approximation is valid for / < kR^ with the further requirement that |/ — kR\ is larger 
than approximately [kUy^^. These conditions have consequences on the accuracy of the 
semiclassical trace formula which is derived from it. The first condition is related to the 
fact that the semiclassical approximation will be given in terms of the classical orbits which 
satisfy / < kR. The consequence of the second condition is that the contribution of or- 
bits to the trace formula becomes less accurate, as the angle of incidence approaches 7r/2 
(whispering gallery orbits). An improvement of the approximation is possible by using a 
uniform approximation for the Hankel functions [23] . 

Inserting the Debye approximation for the Hankel functions and their derivatives into 
eq. (27) we obtain 



Si{k) 



-iifl(k) 



^ K -ipi)e 



-8((7i-7r/4) 



(k + ipi) e+'(<^!-'r/4) 



exp 



2(7; h 2 arctan — 

2 \K 



(29) 



where 



ai = k^R' - P/k^ 



arccos 



— - and 

kR 



Pi 



^k^ - P/R-' 



(30) 



Eq. (29) determines the semiclassical scattering phases (fii{k). We now make use of the 
inside-outside duality in order to obtain from these phases a semiclassical approximation 
for the density of states of the inside problem. From relation (28) it follows that [k > 0) 



oo oo 



d{k) 



/= — oo m = l 



/= — oo 



dcpiik) 



dk 



(31) 



where 6p is a periodic delta- function with period 27r. The wave numbers for the circle 
satisfy / < kR^ and since in this region the derivatives d(fii{k) / dk are positive, the absolute 
value can be omitted from eq.(31). 

Now the delta-function 6p is expressed as a sum over exponential functions. Neglecting 
the constant term in this sum, which contributes to the smooth part of the level density, 
we obtain 



n^pU i = — oo 



n=l 



(32) 



This is equal to the general result for dosc{k) in the scattering approach [15]. 

The trace of the derivative of the n-th power of the S'-matrix is semiclassically calculated 
by applying the Poisson resummation formula 



-I n 



E 



m = — oo 



dLpiik) 



OO dk 



-in(^l {k) — 27Timl 



(33) 



and evaluating the integral over / by a stationary phase approximation. 
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The derivative with respect to / of the arctan-term in eq. (29) is of order 1/kR smaller 
than dcr;/d/ when |/ - kR\ > {kRfl^. It remains smaller in comparison with d(T;/d/ when 
|/ — kR\ approaches [kRY^^ ^ and can therefore be neglected in the region of validity of the 
Debye approximation. The stationary phase condition for / is then given by 



, 7rm 

l-n m = kRsign[m) cos 

n 



(34) 



This condition leads to a restriction of the sum over m to values for which |m| < n/2. 
For |m| = n/2 the stationary point is at / = where d(fii/dl is discontinuous, and the 
integration therefore is evaluated only on one side of the stationary point, which leads to 
an additional factor of 1/2. 

The numbers n and m can be given a physical interpretation. The stationary points of 
eq. (34) correspond to the classical periodic orbits of the circle billiard, which are charac- 
terized by the number of reflections on the boundary n and the winding number m. The 
saddle point ln,m is the angular momentum of the orbit. 

Performing the stationary phase integration and inserting the result into eq. (32) finally 



gives 



oo [n/2] 

dosc{k) = g 

n=2 m = l 




sm 



n 



cos 



. , 7rm\ 

2nkR sm ( 

n } 



-n 1 \- 2n arctan — sm 

2 4 \ K \ n , 



where 



m = n/2 
m ^ n/2 



(35) 



(36) 



This formula expresses the oscillatory part of the spectral density in terms of a summation 
over all families of periodic orbits, which are denoted by the two integers n and m. Its form 
agrees with the general semiclassical formula for integrable systems of Berry and Tabor 
[6, 7], but it has additional phases due to the boundary conditions. The first term in the 



cosine is the action along the periodic orbit kLn^-rm where L,f 



2ni? sin(7rm/n) is the 



orbits length. The phase — 37rn/2 is due to n refiections on the boundary and n conjugate 
points along the orbit. In comparison to the result for Dirichlet boundary conditions there 
is an additional phase that consists of a contribution of 



( k . f 7rm\ 

2 arctan — sm 

\K \ n J 



(37) 



for every refiection on the boundary, where sin(7rm/n) is identical to the cosine of the angle 
of incidence. This is the phase that is given in eq. (9) in the introduction. 



2.2 The Rectangle Billiard 

Quantization of the rectangle billiard with mixed boundary conditions has a twofold pur- 
pose. First, it serves as a check for the generality of the results that were obtained for 
the circle and the semi-circle billiards with mixed boundary conditions. Second, as the 
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derivation unfolds, we arrive at a useful expression that can be used to calculate the semi- 
classical spectral density for a class of higher-dimensional billiards with mixed boundary 
conditions. 

For simplicity, we quantize a X Ly rectangle, with Dirichlet boundary conditions 
on the edges of length and mixed boundary conditions on the edges of length Ly (see 
figure 2). Due to separability, we readily get the quantization conditions 



Dirichlet 




Figure 2: General definitions for the rectangle billiard. 



kj^ YiLx ~1~ 2 arctan 



±1,±2,... 



(38) 



and 



±1,±2,..., (39) 

where k"^ = k1-\- ky and (ibn, ±m) represent the same quantum state. The left hand side of 
(38) is monotonically increasing in k^^n from — oo to +00 and thus guaranteeing a unique 
solution for every n (which justifies the notation). Similar considerations also apply to 
eq. (39). 

The spectral density for the rectangle can be written as 



n,m = l 

00 00 00 

J2 8{E-E^^)- J2S{E-E.,o)- J2 ^(E - Eo-ra) + S{E - E 



n,m = — 00 



m = — 00 



(40) 



where E = k^ , E, 
relations A;^ _„ = —k 



kl ^ + ky ^ and we took advantage of the antisymmetry 



-ky^m which are easily derived from the quantization 
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conditions (38, 39). In particular, £"00 = 0. Apply Poisson summation to the first term 
in (40), using the natural continuations of (38, 39) to real m,n: 



,2 

n.ra = — 00 n.ra = — oo 



00 n 

/ An Am 8 \E — k'^{n) — A;^(m) e 



2^i[pn-\-qm) 



V r Ak, Aky^^S [E-kl- kl) e2-b^fc.)+.™(fc.)] 



p,g=_oo-'-°° Ak^Aky 



1 ^ /'27r 

= - V / d^n'(A;cos^)m'(A;sin^)e2-M'=™«'')+''™('=-"'')] . (41) 

The first change of integration variables (m, n) — )• (A;^;, A;^^) is allowed due to the monotonic- 
ity of n[kx)^ rn{ky) mentioned above. The second change of variables is just transforming 
to polar coordinates {k,6), such that the k integration can be explicitly performed, elimi- 
nating the 8 term. The primes in the above expression denote differentiation with respect 
to the argument. Similar considerations applied to the other terms of (40) finally lead to 

1 ^ /'27r 

S^^\E) = - Y / d^ 4^^)(A;cos^)4]f)(A;sin^)e2-W'=™«'')+''™('=-"'')] , (42) 

8 p,q = -co 

where 

&{k) = n'{k) - 6,,oS{k) , H^) = m'{k) - 8,^,8{h) . (43) 

For p = q = eq. (42) gives the smooth two-dimensional spectral density 

P^\E) = - r A0 S^^\k COS e)Sj^\k sine) , (44) 
S Jo 

where d^^^\h) = n'[k) — 6[k) and d^y^\h) = m'{h) — 6[k) are the one-dimensional smooth 
spectral densities (the (5-functions are a convenient way to represent the constant term 
of N^^^^kj). The "combination formulae" (42, 44) thus express the two-dimensional 
spectral density using the one-dimensional components in a simple way and rely on the 
energy decomposition relation E = k^ -\- ky and the positive-negative antisymmetry of 
kx^±m ky^±m- Consequently, it can be extended to calculations of the spectral density for 
higher-dimensional billiards, provided that the energy can be decomposed in the above 
way. In particular, it is useful for three-dimensional "cylindrical" billiards [21], i.e. , two- 
dimensional billiards elongated along the perpendicular axis and closed from above and 
below. For example, the cubic billiard with a variety of mixed boundary conditions can be 
easily calculated this way, including the smooth contributions of right angle edges and cor- 
ners, which are otherwise difficult to estimate using more general methods. If one assumes 
that the Weyl series has the general form with volume term, surface term et cetera, also in 
the case of mixed boundary conditions, then these results can be applied also to billiards 
other than the simple ones mentioned above. 

To derive semiclassical expression for S^^\E)^ we use the saddle point approximation 
to evaluate the oscillatory integrals {{p^q) 7^ (0,0)) in (41). The phase appearing in these 
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integrals can be divided into two parts: 



^P^e) = 2Tr [pn{k cos 9) + qm{k sin 9)] 
= 2 ^ ^ 

k cos 9 



( k cos 9\ , 

pkLx cos 9 -\-2p arctan + qkLy sm 9 



ip^^{9) + 4p arctan 



(45) 



where (^^^{9) is the phase that appears in the pure Dirichlet case and is rapidly oscillating 
in 9 for large k. The phase (f>^{9)^ which is genuine to the mixed case, is both bounded and 
slowly oscillating and thus can be appended to the slowly varying prefactors as discussed 
in subsection 2.1 (see eq. (34)). Thus, the saddle points are the same as for the pure 
Dirichlet case: 

ii,2 _ P^y 



tan 9 



pq 



(46) 



The angles 9^'^ correspond to the angles between the x-axis and classical periodic or- 
bits (tori) which trace the x dimension p times and the y dimension q times (see fig- 
ure 2). Putting together the above results, we finally get the semiclassical approximation 
to d(^^\E): 



J(2^)(i?) = J(2^)(i?) + 42Z))(^^ 
La: Lt, 



LxLy 
47r 



^irk 4:Trk 



+ 



L, 



2TTk\ l + 



+ 



LxLy 

47r 



E 



Lx 
2Trk 



CO 

I] cos 

p=l 



TrkL 



cos 



pq 



27rF j^l + 
'k 



kLpq -\- 4:\p\ arctan ( — cos 9pq 



2kpLx + 4p arctan 



L 

^j:cos{2kqLy 

q = l 



(47) 



pq 



where Lpg = 2^{pLxY + [iLyY is the length of the [p^q] periodic orbit and cos ^, 
I cos^p^^l = 2\p\Lx/ Lpq. In the smooth part we kept terms of all orders, since we evaluate 
the integrals exactly, but in the oscillatory part we kept only leading order terms to be 
consistent with the saddle point approximation. The double sum relates to the periodic 
tori, while the two single sums relate to "boundary orbits" [9] that bounce along the edges 
of the rectangle. Examining the above expression, we conclude that the effect of the mixed 
boundary conditions on the oscillatory part is simply to modify the contribution of each 
periodic orbit by a phase, which is 2 arctan [{k/n,) cos9] for each bounce from a wall with 
mixed boundary conditions, where 9 is the angle of incidence with respect to the inside- 
pointing normal. The same rule also applies to the boundary orbits, for which 9 is either 
or 7r/2. This result, as well as the modifications to the smooth spectral density (length 
term and right angle Dirichlet-mixed corners) are equivalent to and consistent with the 
results that were obtained above for the circle (see (24,37) ). 
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3 General billiard systems 



We now consider general smooth two-dimensional billiards. A very convenient method 
for the treatment of billiard systems are boundary integral methods since they reduce 
the two-dimensional billiard problem to an effectively one-dimensional problem along the 
boundary of the billiard. However, as has been discussed in the introduction, in the case 
of mixed boundary conditions one encounters the difficulty that the integral kernels have 
a singularity, which makes a direct perturbation expansion impossible. One possibility to 
solve this problem, is to apply a transformation to the integral equation, which leads to an 
non-singular kernel. This corresponds to a resummation of the equation in all orders of n 
and we will use this method in order to derive the smooth part of the level density [3]. For 
the oscillatory part of the level density, however, we choose two different approaches that 
avoid the problem of singular integral kernels. We first apply methods of scattering theory 
in order to obtain the semiclassical approximation for convex chaotic billiard systems. For 
the calculation of the scattering matrix the Kirchhoff approximation (presented in [15]) is 
extended to mixed boundary conditions. We further derive the semiclassical approximation 
for a concave billiard, the Sinai billiard, by the KKR-method. 

3.1 The smooth part of the spectral density 

In reference [3] Balian and Bloch derive the asymptotic form of N[E) for an arbitrary 
smooth three-dimensional billiard with mixed boundary conditions. We use their method 
in order to obtain the corresponding two-dimensional result. Our notation differs from 
that of Balian and Bloch in that the Green function and the normal derivative have a 
different sign. 

Consider a two-dimensional domain D with area A and a smooth boundary E of length 
L. The Green function of the Helmholtz equation for this domain is determined by the 
differential equation 

{V'' + E)G{f,f',E) = 6{f-f') (48) 

with the boundary condition 

{K + dn)G{f,f',E) = 0, rGS. (49) 

It is related to the smoothed spectral density d^[E)^ which has Lorentzian peaks of width 
p at the energy values by 

dJE) = -i / dV [Im Gif, r ', E + z//)]//=/ . (50) 

TT J A 

For the derivation of the large A;-behavior of the Green function the differential equation 
(48) and boundary condition (49) are replaced by an integral equation. This is done by 
representing the Green function G'(r, r*', E) as the sum of the free Green function and a 
single layer potential: 

G{f, r ') = G'o(r, r ') + ^ds„ G'o(r, fc)ji{fc, r ') (51) 
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with a density fi[fp^f') which is determined by 

IJ.{ff3, r ') = 2{k + dn^) Go{ff3, r ') + 2 /^ds„ [k + dn^]Go{ff3, fa)iJ.{fa, r ') . (52) 



s 



For brevity of notation, the energy dependence of the Green functions and the density fi 
will be omitted, and Greek indices will denote coordinates on the boundary. 

For Neumann boundary conditions [k, = 0) the integral equation (52) can be solved by 
iteratively replacing the potential fi in the integral term of eq. (52) by the whole expression 
for /i. For k, ^ this cannot be done since the Green function Go{fp^ r*„), in contrast to its 
normal derivative, has a (logarithmic) singularity at r„ = fp. For that reason, the integral 
equation has to be transformed into a different integral equation, for which the kernel 
is uniformly bounded. This is done by introducing an auxiliary Green function r(r„,r*^) 
which is defined on the boundary E by 



r(r„, ffs) - 2k y^ds^ G'o(r„, r^) r(rl,, ffs) = 6{fa - rp) . (53) 
Multiplying eq. (52) by r(r„,r*^) and integrating over rp one obtains 

+2 j dsp ds^ r(r„, fp) dn^Go{fp, f^) ji{f^, r ') . (54) 

This equation is now solved by a perturbation expansion and the result is inserted into 
eq. (51). One obtains 

G{f, r ') = Go{f, r ') + 2 J^ds^ dsp Go{f, f^) T{f^, fp) [k + 9^,^) Go{fp, r ') (55) 

+4: j^dsadspds^fdssGo{f,fa) T{fa,fp) dn^Go{fp,f^f) T{f^f,fs) [k + dn^) Go{fs,f') + . . . 

Inserting this into eq. (50) yields the following expression for the level density: 
A 1 If 

d^{E) = Im 2 dsc,dspV{fc,,fp){n + dng)F{fp,fc,) (56) 

47r TT L 



+ 4 y^ds„ dsp ds^ dss r(r„, r^j) dn^Go{fp, f^) V{r^, n) [k + dn^) F{fs, r„) + . . . 
where the function F[fa^fp) is defined by 

F{ra,fp)= / dVG'o(r„,r ) G'o(r, r/j) , (57) 



and the energy of the Green functions is given hy F -\- iv . In contrast to the three- 
dimensional calculation, the introduction of a convergence factor is not necessary here 
since the integrals are convergent. 
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The perimeter term 



The evaluation ol the integrals in eq. (56) is based on the fact that the functions Go and F 
are short-range functions when the energy E is large. In a first approximation it is assumed 
that the range is smaller than the smallest radius of curvature of the boundary E and the 
boundary is locally replaced by its tangent. In this approximation only the first integral 
term in eq. (56) contributes, since dni^Go{fp^f^) vanishes if the curvature of the boundary 
is zero. One proceeds in the following way: the coordinates r„ and fp are replaced by 



rf3 



(58) 



and the integral over f(_ is evaluated along the tangent of the boundary at the point 
Along this tangential line the functions Go and F are translational invariant and for that 
reason it is convenient to introduce the Fourier transform along the tangent: 



fip) 



Ztt 



(59) 



In terms of the Fourier transforms the contribution to the level density can be written as 

A 2 



dAE) = - ^Im 

47r TT^ 



ds„ 



Us 



ApV{p) (KF{p) + dnF{p) 



+ ... 



(60) 



The calculation of the Fourier transforms of the functions Go, F and F is done exactly as 
in the three-dimensional case of Balian and Bloch. For that reason, we give here only the 
results: 

1 A. X a[p) 



Go{p) 



2a{p) 



T{p) 



K + a{p) 



(61) 



where a[p) = \Jp^ — , ki, = \/F -\- i 



V anc 



) 27r a?{j)) + 



Fip) 



8a^{p) 



dnF{p) 



Sa'^{p) 



(62) 



Those of the functions in the above equations which depend only on one parameter are 
the Fourier transforms of functions with arguments in the plane, which are translationally 
invariant along the plane. Go{p,z) is the Fourier transform of Go(r„,r*) and z is the 
distance of r from the plane. The expressions in equations (61) and (62) have exactly the 
same form as in the three-dimensional case, the difference being that p is now a scalar 
instead of a two-dimensional vector. The results (61) and (62) are inserted into eq. (60) 
and the limit ;/ — )• is performed, yielding the first asymptotic terms for the mean level 
density. Multiplying the result by 2k gives the result for the mean level density in terms 
of the wave number k: 



d{k) 



A , Lk 
-k 



2tt 27r2 



2tt Att 



1 



K — a 



Im / da 



+ ... 



+ 



l + (f 



(63) 
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The curvature term 



The next term in the asymptotic expansion of d[k) is obtained by including those cor- 
rections to the tangent approximation which are linear in l/i?c^, where R^^ is the radius 
of curvature at r^^. Such corrections exist for F[fa^fp)^ dn^F[fa^fp) and dn^Go{fa^fp)^ 
whereas for Go{fa^fp) and r(r„,r*^) the first corrections are quadratic in 1/Rco and can be 
neglected. The integrals in eq. (56) are expressed in terms of the Fourier transforms of the 
Green functions and their corrections and the next asymptotic term for the energy level 
density follows as 



— rim 
— rim 



ds^ / dpT{p)[KSF{p) + SdnF{p)] 
Ue Jo 

oo ^ ^ ^ ^ ^ 

d5,, / dpT{p)6dnGo{p)T{p)[KF{p) + dnF{p)] 







(64) 



The Fourier transforms of the Green functions and their corrections can again be derived 
exactly as in the three-dimensional case [3] and give the result 



SdnGoip) 



6F{p) 



a'^jp) + kl 
' m^a^{p) 



SdnFip) 



2a^ip) + Ul 
IGR^a^p) 



(65) 



These expressions are slightly different from their counterparts in three dimensions, since 
their derivation involved the divergence of p which depends on the dimension in which it 
is calculated. The terms in (65) are inserted into (64), the limit ;/ — )• is performed and 
the result is multiplied by 2k in order to obtain the next term in the asymptotic expansion 
of d[k). The result is 



AT 
d(k) = —k 



Im 



27r 47r 



di 



da 



a^K^ + klK^ -a^- 2a^kl 



eR^ J-^K a^{a + kY J a"^ + kl 

2 1 
1 — , + 



i + (f 



+ ... 



l+(f 



iirk 



l + (f 



-1 3/2 



3/2 



di 



(66) 



+ ... . 



An integration over k gives 



N{k) 



47r 47r 



'1 + 



127r 



di 



+ ... 



(67) 

This is the final result for the leading three asymptotic terms, the area, perimeter and 
curvature term, for the mean spectral staircase for an arbitrary smooth two-dimensional 
billiard with mixed boundary conditions. For the determination of the integration constant 
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in the integration of eq. (66) see the dicussion in the section on the circle biUiard. Eq. (67) 
is identical to eq. (25) when applied to a circle. Furthermore, in the limits /€ — )• oo and 
/€ — )• eq. (67) gives the correct expressions for the case of Dirichlet and Neumann boundary 
conditions. 

3.2 The oscillatory part of the spectral density 

The oscillatory part of the spectral density for chaotic billiards with mixed boundary 
conditions (4) is considered in this subsection. We use the scattering approach, in which 
the billiard is viewed as an obstacle for a scattering problem and use (32) for dosc{k). The 
scattering matrix is obtained using the Kirchhoff approximation in a way which applies to 
strictly convex billiards. The approximation assumes that the incoming wave locally sees a 
straight line at the scattering point. The case of Dirichlet boundary conditions is given in 
detail in [15]. For mixed boundary conditions the additional phase (9) is again recovered. 
This phase is slowly varying and in saddle point approximations its value is therefore taken 
at the saddle point (in section 2.1 this point is justified for the circle). 

We consider the scattering solution which consists of an incoming plane wave V'm(^) = 
gjfcj-r (-with momentum ki and direction 9i) and a scattered wave il'scati'f^)- The scattering 
amplitude /(^i. Of) to an outgoing direction ^/ is connected to the asymptotic form of the 
scattered wave at large r by 



The momentum in the outgoing direction is denoted by kj^ so that k = \ki\ = \kf\. Using 
Greens theorem the scattering amplitude is expressed as a boundary integral: 



This expression is exact, but involves the knowledge of both tf^scatif^) and dnt/jscatif^) on the 
boundary. The boundary condition (4) gives one linear combination of these two quantities, 
another is found by using a short wave-length approximation. 

Given an incoming direction the billiards boundary is divided into an illuminated 
part (Eg) and a shaded part (S^), dividing the expression for the scattering amplitude (69) 
into two parts, / = fi + fs- For the shaded part the approximation is introduced by arguing 
that the wave function is vanishingly small close to and hence t/j[f) = dnt/j[f) = 0. 
This results in exactly the same expression as for Dirichlet boundary conditions, giving a 
semiclassical approximation for the forward scattering peak [15]. The scattering matrix (in 
the angle representation) is then connected with the illuminated part only. 

The short wave-length approximation assumes that the behavior of the scattered wave 
at a point on the illuminated boundary is determined by specular refiection off that point 
and that the surface is locally a straight line. 

We first consider a wave scattering off a straight wall, with a normal h pointing towards 
the side where the scattering process is taking place. The incoming wave is given, as before, 
by V'm(^) = e^^''^ . The scattered wave at any point r*on the wall and its normal derivative 





ds \dn4)scat{r ) + 1 (kf ■ h{f) ) tpscat{r)\ e 



(69) 
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at that point, are expressed in terms of the incoming wave by 



i^scatir ) = -e-'^tP.^f) , (70) 

dn^hcaiir ) = i{h ■ n)e-''^^ha{r ) . (71) 
Imposing the boundary conditions at the wall, the additional phase (f> is determined to be 

^ = 2 arctan ( — cos ^ ) = 2 arctan I — -| , (72) 

\K J V " / 

where 6 is the angle of incidence and is just the additional phase (9) which appears for 
mixed boundary conditions. Note that (71) is valid for a wall at any angle or position. 

We use the result for the straight wall in (69) to obtain the approximation for the 
scattering matrix 

S{0^, 0f) = :^M0^, Of) = ^ [ ds {% - kj) ■ n{fyi^'-hyr-^m , (73) 
Ztt 47r Js, 



where 



/ ki ■ h(f ) \ 

(r) = 2arctan ^ . (74) 



The integral (73) is evaluated using a saddle point approximation, where the additional 
phase (f>{f) is taken at the saddle point. The saddle point Tq for a general convex billiard, 
is the point of specular reflection for a trajectory with incoming direction 6i and outgoing 
direction Oj. The scattering matrix is finally given by 

S{0.. 0f ) = -l (^) ' [{kf - ki) ■ n(ro)i?(ro)] ' e'-(^.-^,)-^^o-.^(/o) ^ (75) 

where R[f) is the radius of curvature and (f>{f) is given by (74). Note that, since r*o 
is the point of specular refiection, we have kj ■ n[fo) = —ki ■ n[fo) = A;sin(|^j — 6i\/2). 
The scattering matrix for the circle may be directly obtained using R{fo) = R and Tq = 
Rh[fo). By transforming to the angular momentum representation, using a saddle point 
approximation, we recover (29). 

The oscillating part of the spectral density, given in the scattering approach by (32), is 
evaluated using the saddle point method, and the additional phase (f> is taken at the saddle 
point. The calculations for Dirichlet boundary conditions [15] may then be followed, with 
the appearance of the additional phases being the only difference. Assuming that all 
periodic orbits are isolated, we obtain the Gutzwiller sum 

1 ^ ^ cos \m(kL^ — z/^f- — n^TT + $^) 

doscik) = - E E L, ^ ; ' ; , (76) 

^ 7 i^i ,/2-TrM™ 



where the sum is over primitive periodic orbits 7 with repetitions m. Each primitive orbit 
has a length and a monodromy matrix M^. The number of refiections along the orbit is 
and is the maximal number of conjugate points along the orbit (the Maslov index). 
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The only difference from tlie formula for Diriclilet boundary conditions is the appearance 
of the additional phase which is just the sum of the additional phases accumulated 
along the orbit: 

$^ = ^(^; = 2^arctan Pcos^;] , (77) 

i i \ J 

where 6^^ are the angles of incidence at the reflection points along the orbit which are 
labeled by i. 



3.3 The Spectral Density for the Sinai Billiard 

In this subsection we describe very briefly the semiclassical quantization of the two- 
dimensional Sinai billiard (disc of radius R embedded inside a square of side L) with 
mixed boundary conditions (4) on the disc. For this purpose we follow Berry [8] and use 
the KKR method [24, 25, 26]. We first unfold the billiard into an infinite lattice and con- 
sider a general potential with circular symmetry and range R in the center of each cell 
(2i? < L). Using the Green theorem and taking advantage of the translational and circular 
symmetries, we arrive at a set of linear equations, that have a non-trivial solution (i. e. , 
an eigenvalue of the billiard) if 



({h) = det 



6u, - tan ( ^^^4~^ I (^u'ikL) 



, (78) 



where a;;/ are "structure constants" of the lattice that are independent on the potential and 
^i{k) are the scattering phase shifts that include all the information about the potential. 
For mixed boundary conditions the phase shifts are given by 

S = = _ '^Hr{kR)-kHr'{kR) 
' KH+{kR) - kH+'{kR) ' ^ ^ 

which is different from eq. (27) for the circle, due to dn = —dr in the Sinai billiard case. 
Semiclassically we have 

^i{k) = ^^{k)-Mk). (80) 

where '~pf{k) are the phase shifts for Dirichlet boundary conditions on the disc and = 
2 arctan [{k/ k.) cos 6i]. 6i is the angle of incidence to the disc of a classical particle having 
angular momentum /. 

To calculate the expression for the spectral density d{k) one considers ( — l/vr) log C(A;) 
{C{k) can be made real [8]). Then one uses the matrix relation log det = Tr log and 
formally expands the logarithm in a Taylor series. To get the semiclassical approximation 
of d{k), Poisson summation is used and the resulting integrals are evaluated by using saddle 
point approximation term by term. (The detailed and rather long derivation is one of the 
main subjects of [8].) Substituting in (78) and performing these steps, we immediately 
get a semiclassical expression for the Sinai billiard with mixed boundary conditions that 
contains no expansions in powers of n. Thus "/€ is resummed to all orders" by the KKR 
determinant and is included only through the phase shift cf)i. Moreover, due to the slow 
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variation of (pi as explained above (see subsection 2.1), the modifications are particularly 
simple and result in 



k{L'' - 7ri?2) L R 



27r 



TT 



+ dbb{k) 



cos 



7 m = l 



mkL^ — mn^ 



TT + m ^ 2 arctan [ — cos 91^ 



8 = 1 



f8f^ 



where A-^^m are the semiclassical amplitudes that appear in (76) and n-^ are the number 
of collisions of the primitive periodic orbit 7 with the billiard walls. Note the fact that 
for the Sinai billiard the Maslov index z/^ is always 0. The inner sum is over collisions 
with the disc [9^^ are angles of incidence) and contains the modifications to dosc{E) due 
to mixed boundary conditions. The ^-independent term dbb{k) contains the non-generic 
contributions of the "bouncing ball" orbits [8, 9, 27] and will be further discussed below 
(section 4). It is interesting to note, that although in eq. (80) the sign of the additional 
phase shift (pi is opposite to the that of the circle (and more generally the convex) billiard, 
the final results for both d[E) and dosc{E) agree in signs with the previous cases, which 
emphasizes the generality of the modifications implied in eq. (81). 



4 Application of Mixed Boundary Conditions to the 
Spectral Diagnostics of the Sinai Billiard 

The spectral density of the Sinai billiard, in either two or three dimensions, and in the 
presence of mixed boundary conditions on the inscribed scatterer (either disc or sphere), 
can be written as 

d{k] 6, a) = d[k; 6, a) + dosc{k] 6, a) , (82) 

dosc{k; 6, a) = dgen{k; 6, a) + dbb{k) + dtan{k; 6, a) (83) 

where d is the smooth density, dgen comes from contributions of generic, isolated and 
unstable periodic orbits (see (81)), dbb represents the contributions of the non-generic 
"bouncing ball" (BB) families of marginally stable orbits that do not collide with the 
scatterer [9, 8, 27], and dfan contains the non-generic contributions due to periodic orbits 
which are tangential to the scatterer (see figure 3). When one is interested in the generic 
part of the spectral density, one must find a reliable and accurate method for the elimination 
of the BB and tangent contributions. This is not an easy task, since the contribution of the 
leading-order BB families are 0[k^^^) stronger than that of an isolated periodic orbit in the 
two-dimensional case [9], and 0{k^) stronger in three dimensions. In the semiclassical limit 
this may amount to a very large factor, so that to identify the contributions of a single, 
generic orbit one must subtract the BB contributions in an accurate and reliable way. This 
problem is especially acute when one considers the cosine transform of the oscillating part 
of the spectral density - the length spectrum: 

D{x;k^ax;b,a) = / dk dosc{k;b,a) cos{kx) . (84) 
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Here, due to the finite range of k values, tlie lengtli spectrum is composed of peaks witli a 
finite widtli, and tlie contribution of tlie BB families might overlap and obscure the peaks 
due to unstable, isolated periodic orbits. This problem is very well demonstrated in the 
top part of figure 4. At the length x = 2 there exist a BB family and an isolated periodic 
orbit (a two fold repetition of the shortest orbit of length 1). The latter has an amplitude 
which is smaller even than that of the peak at x = 1. 

In order to eliminate the non-generic contributions, it is natural to seek analytic ex- 
pressions, and to subtract them from the spectral density. Indeed, for the two-dimensional 
stadium billiard the contribution of the (single) BB family were calculated by Sieber et 
al. [9], and were successfully applied in the spectral analysis. These calculations become, 
however, much more intricate for the Sinai billiard, in particular for three dimensions. In- 
stead of a single BB family as in the stadium, the Sinai billiard admits a rich variety of 
BB families, especially in three dimensions [28]. Moreover, in three dimensions the total 
(geometrical) measure of each family of BB is difficult to evaluate. Tangent orbits (which 
are not encountered in the stadium billiard) require a special treatment, which is not yet 
available. 

Facing all these difficulties, we were forced to introduce a new method for the elimination 
of the BB and tangent orbits. It makes use of the extra freedom which is gained by the 
variation of the boundary conditions. As a matter of fact our interest in the semiclassical 
theory of billiards with mixed boundary conditions stemmed originally from this particular 
application. 

The main idea of this method is to look at differences of spectral densities obtained 
with boundary conditions which differs only on the scatterer. Since the BB orbits do not 
collide with the scatterer, they are insensitive to changes of the boundary conditions on 
it. Therefore, their contributions will be eliminated in the difference. We shall discuss 
here two convenient schemes. In the first, we shall subtract the spectrum with Dirichlet 
boundary condition on the scatterer from the spectrum with the Neumann condition on 
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the scatterer. The difference wiii be caiied the N-D spectrum. 

Talcing into account the explicit expressions for dose for Neumann and Dirichlet bound- 
ary conditions we get from eq. (81) 




(85) 



where is the number of collisions with the straight boundaries, and the prime over the 
sum indicates that only periodic orbits with an odd number of collisions with the scatterer 
should be considered. d[^^ is the residual contribution of the tangent orbits. The BB 
contribution dbb[k) is completely eliminated by the subtraction. This result is valid in 
the semiclassical approximation, and its applicability in practice should be tested. This 
was done by analyzing numerically the two-dimensional quarter of a Sinai billiard with 
Z = 2, i? = 0.5. Our numerical database consists of the lowest 5700 energy levels [29] in 
the A;-range < kn < 300, for both Dirichlet and Neumann boundary conditions on the 
disc. 

In figure 4 (top) we present the two length spectra for Neumann and Dirichlet boundary 
conditions on the disc and one clearly observes the large contributions of the BB families, 
which have lengths of 2 and 2\/2 [8]. The change of sign between Dirichlet and Neumann 
boundary conditions due to odd number of bounces on the disc is evident, e.g. for the 
shortest orbit of length 1 which bounces only once from the disc. There are also peaks 
that do not change sign, as expected. The N-D spectrum is shown in figure 4 (lower part). 
The peaks corresponding to the BB are greatly diminished, and generic periodic orbits can 
be distinguished more easily. There are, however, some drawbacks to the N-D analysis. 
The first one is that it eliminates half of the generic contributions to the spectral density, 
namely those orbits that bounce an even number of times from the scatterer (see eq. (85) ). 
This can be easily rectified by comparing spectra with other boundary conditions on the 
scatterer rather than the special N-D analysis. A more serious drawback is due to the 
contributions of tangent orbits. They have the same length as the BB, and their residual 
contribution can be clearly seen in figure 4 for x = 2, 2^2. It was found numerically to be of 
the same order in k as for isolated orbits. This gives an indication, that in three dimensions 
such tangent orbits, that form 1-parameter families, are likely to give contributions of order 
A;^/^, which can obscure a significant part of the generic contributions. 

The second, and more powerful comparison method consists of taking the derivative of 
d[k; 6, a) with respect to a at a = 0. Using eqs. (9, 81, 83) we get 



The term dbj, dropped trivially, and the oscillatory contributions are now multiplied by the 
new prefactors J2iZi cos 91^. These prefactors vanish (only) for tangent orbits {91^ = 7r/2, i = 
1, . . . , n'i). If we now assume that the tangent orbit contributions behave similarly to those 
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Figure 4: Top: length spectra of quarter of the two-dimensional Sinai billiard, Z = 2, i? = 
0.5, with Neumann and Dirichlet boundary conditions on the disc. Bottom: difference 
between Neumann and Dirichlet length spectra. Isolated periodic orbits are indicated by 
vertical lines and bouncing ball manifolds are indicated by daggers. 
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of the periodic orbits, it will follow that their contributions will also be eliminated. Hence, 
the a derivative of the spectral density consists of generic contributions, exclusively. 

To examine (86) we use instead of the plain length spectrum (84) a Gaussian modified 
one with the kernel 

j{k, X] ko, a) = J—^ exp (- ^^ ) cos{kx) . (87) 

The parameters ko and a were chosen, such that for the calculated spectra, the integration 
limits can be safely taken to iboo, thus avoiding undesired boundary effects. Integrating 
both sides of (86) with this kernel, one gets 
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Here we used the explicit expression for the smooth spectral density of the two-dimensional 
quarter Sinai billiard including the corner and curvature terms. To check eq. (88), we 
calculated numerically the derivatives of energy levels with respect to a by using finite 
differences for a = and a = 0.0005 7r/2. In both cases b = 200 and < kn < 300. 
In figure 5 we present the quantum results together with the semiclassical predictions of 
eq. (86) which consists of contributions due to generic isolated and unstable periodic orbits 
(and their repetitions) only. The agreement between the quantum and the semiclassical 
results is good throughout the entire length interval (but for a few deviations which are 
due to special effects which will be discussed latter). This demonstrates the power of the 
derivative method, and the general success of the semiclassical theory. 

A more detailed examination of figure 5 reveals the following details. At the lengths 2 
and 2\/2 which correspond to the shortest BB families and their limiting tangent orbits, 
we can observe very small structures. Comparing them to the corresponding peaks in the 
N-D spectrum, we see that the tangent orbits behave as expected. The very small peaks 
that nevertheless exist there, are attributed to higher-order corrections to the semiclassi- 
cal approximation and to the small contribution of the second repetition of the shortest 
periodic orbit of length I. Other BB families cannot be separately examined for the given 
k range due to the clustering of periodic orbits (the total number of periodic orbits grows 
exponentially with their length). 

It is important to note that the quality of the semiclassical theory is uniform through- 
out the length domain, even at longer lengths, where individual periodic orbits cannot be 
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resolved. This is a very demanding test for the semiclassical theory since it involves the in- 
terference of dozens of periodic orbits per wavelength. To get such a good reconstruction of 
the exact length spectrum requires a very accurate evaluation of the individual amplitudes 
and phases. 

The largest discrepancy between the numerical data and the semiclassical theory is 
seen in the vicinity of x = 6.21. At this length there exists a triplet of unstable periodic 
orbits which are very close to tangency with the disc. It would be eliminated by pruning 
if the radius of the disc was increased beyond R = 0.555. Other points in the length 
spectrum where deviations of the semiclassical theory from the numerical length spectrum 
can be observed are also close to bifurcations of this type with two or three periodic orbits 
becoming simultaneously tangent to the disc. Therefore we suspect that the deviations 
have their origin in the failure of the semiclassical theory in its present form in the vicinity 
of tangent orbits. 

The method presented above was recently applied in the analysis of the 3-D Sinai 
billiard. There, the length spectrum is entirely overwhelmed by BB contributions. One 
can study the contributions of isolated periodic orbits only after the spectrum is analyzed 
by the methods explained above. The detailed analysis of the 3-D Sinai billiard will be 
presented elsewhere. Finally, we would like to note that the method presented here can 
be applied in many variations, whenever a particular set of orbits is of particular interest, 
and one wishes to suppress all other contributions. 

5 Parametric Spectral Statistics 

Random matrix theory has been successful in describing universal properties in the spectra 
of a large variety of single systems. These include classically chaotic systems as well as 
disordered mesoscopic systems. There is now an increased interest in considering not 
only single isolated systems but also systems that depend on an external parameter and 
in describing correlation functions in this parameter and statistics that depend on the 
sensitivity of the energy levels to the change of the parameter [10, 11, 12, 13]. 

Two very often considered quantities are the velocity and the curvature distributions, 
i. e. the distributions of the first and second derivatives of the levels with respect to the 
parameter. For these quantities several theoretical results exist. The velocity distribu- 
tion P{v) is a Gaussian for random matrix ensembles. It has been shown by Gaspard et 
al. [10] that the asymptotic distribution for large curvatures is universal. Subsequently, this 
behavior was verified numerically for several systems [30, 31, 12]. Formulae for the com- 
plete curvature distributions for the GUF, GOF and GSF of random matrix theory were 
conjectured by Zakrzewski and Delande [12] based on detailed numerical examinations. 
In addition, it was suggested that deviations from these distributions in chaotic systems 
could be taken as a measure of the degree of scarring in this system [31, 12]. It has been 
proved recently that the conjectured distributions are indeed the exact distributions for 
the random matrix ensembles [32, 33]. 

In this section we examine the velocity and curvature distributions for the Sinai billiard 
with mixed boundary conditions on the dispersing arc. As has been mentioned in the 
introduction, it is a special feature of this parameter dependence that it restricts only to 
the quantum system, whereas the corresponding classical system is parameter-independent. 
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Figure 5: The derivative of the length spectrum of the quarter Sinai biUiard, Z = 2, i? = 0.5, 
with respect to a at a = 0. Lengths of periodic orbits are indicated by vertical lines, lengths 
of BB families are indicated by daggers. 
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As parameter we take the quantity a of eq. (5) with fixed h = 200. The velocities are 

^"^^ , m = (('^X) , (90) 



d 



a 



where Xn{a) = N[kn{a); 6, a) are the unfolded energy levels. We scaled the first derivatives 
dxn/da by their local standard deviation, since this standard deviation is not constant 
with increasing energy. The brackets < • >x denote an average over the levels Xn in a local 
interval around the value x. The scaled curvatures are defined as 

where p classifies the corresponding random matrix ensemble. It is equal to one for GOE, 
which is appropriate for systems with time reversal symmetry. The distribution of curva- 
tures for GOE is given by 

^('' = sittV • 

We compare the numerical distributions for the Sinai billiard to the random matrix 
results and discuss the deviations that occur, especially in the context of the scarring of 
wave functions in the system. The calculations were done for the quarter Sinai billiard 
which has two symmetry classes, odd and even. The eigenfunctions of these symmetry 
classes are solutions of the Helmholtz equation for the fundamental region of the Sinai 
billiard, which is drawn in figure 7 with thick lines, with Dirichlet (odd) or Neumann 
(even) boundary conditions along the line y = x, mixed boundary conditions along the 
circular part of the boundary and Dirichlet boundary conditions along the remaining two 
straight sections of the boundary. As in the previous section i? = 0.5 and L = 2 and for 
these parameter values there exist two families of bouncing ball orbits. Eor the numerical 
analysis the energy levels with a wave number in the range 100 < k < 300 were used. 
This corresponds to approximately 2500 levels in each symmetry class. The statistical 
distributions were evaluated by averaging over the two symmetry classes. 



5.1 Non-Generic Features in the Parameter Dependent Energy 
Spectrum 

Eigure 6 is a plot of the wave numbers kn of a part of the spectrum as a function of the 
parameter a for the two symmetry classes of the quarter Sinai billiard. There is a mean 
decrease of the levels with increasing a which is due to the a-dependence of the mean 
number of energy levels N[k). A striking feature in figure 6 is that one can distinguish 
apparent horizontal lines. These lines consist of almost straight pieces which are interrupted 
by avoided level crossings. As will be argued below these lines can be attributed to the 
existence of the families of bouncing ball orbits. We discuss these lines in more detail since 
they lead to deviations in the velocity and curvature distributions from the random matrix 
theory. Similar apparent lines have also been found in various other systems (see e. g. [34]) 
and also in experiments [35]. 

The lines are regularly spaced and they can be divided into two families. The first 
family appears at the same positions in the spectra of both symmetry classes. The lines 
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Figure 6: The wave numbers kn{a) in the range 190 < kn < 200 as functions of the 
parameter a for the desymmetrized Sinai biUiard, for odd (left figure) and even (right 
figure) symmetry. The triangles indicate eigenvalues of two rectangular billiards (dark 
triangles - first family, empty triangles - second family), as discussed in the text. 
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actually appear as double lines where the upper line is less pronounced than the lower line. 
The positions are roughly given hy k = nvr, where n is a natural number (dark triangles 
in figure 6). In the case of the second family of lines there is a difference between the 
two symmetry classes. One can distinguish horizontal lines at the approximate positions 
k = n7r/-v^, but they appear only for odd integers n in the odd symmetry class and for 
even integers n in the even symmetry class (empty triangles in figure 6). Although we 
did not calculate wave functions, we strongly suspect that the horizontal lines correspond 
to wave functions that are mainly concentrated (scarred) along one of the two families 
of bouncing ball orbits and that these states are not strongly infiuenced by changing the 
boundary conditions on the circular part of the boundary. Under these assumptions the 
actual position of the lines and the difference between the two symmetry classes can be 
explained by a very simple model, namely by approximating the wave functions by the 
superposition of two eigenfunctions of appropriate rectangular billiards. 

We consider first the position of the lines in the vicinity of a = 0. Later on, we shall 
show that the lines are not strictly horizontal, but that they drop slightly as a approaches 
7r/2. We start with the first family of lines. It is convenient to discuss them in the quarter 
Sinai billiard instead of the fundamental region. We assume that the wave functions that 
correspond to these lines are scarred along the first family of bouncing ball orbits. These 
are orbits that are parallel to the x-axis or the y-axis in the coordinate system shown in 
figure 7, i. e. they occupy the region of the upper shaded rectangle in figure 7 and the 
region of another rectangle that is obtained by refiecting the first rectangle on the line 
y = X. We then express the wave functions by a linear superposition of a function /(x, y) 
and its mirror image /(j/, x), where /(x, y) is concentrated along the first of the rectangles 
only. The wave functions of the two symmetry classes are then given by /(x, y) =p /(j/, x). 
This decomposition in terms of f[x^y) and f[y^x) is of course only possible within an 
approximation since it implies that both symmetry classes have the same eigenvalue. 

The fact that the eigenvalues of the wave functions are only little infiuenced by the 
change of the boundary conditions on the circular part of the boundary points to a small 
excitation of the function f[x^y) in x-direction and a high excitation in j/-direction. For 
that reason, we determined its eigenvalue by an adiabatic approximation in which the 
high excitation in j/-direction and the low excitation in x-direction are separated. In the 
energy range considered, the results differ only little from the eigenvalues of the rectangular 
billiard in figure 7 with Dirichlet boundary conditions. For simplicity we therefore restrict 
the discussion to the eigenvalues of this rectangular billiard. The wave numbers of its 
eigenvalues are given by kn^m = vry^n^/a^ + where n and m are integers and a = 1 

and b = 1 — R are the lengths of its sides. The wave numbers kn^m provide a more detailed 
description of the positions of the lines than the rough estimate k = mr. The values of kn^m 
are marked in figure 6 for n = 61 to 63 and m = I and 2 by black triangles. As one can 
see, they agree well with the position of the horizontal lines. Scarred wave functions that 
correspond to higher values of m might also exist, but they are more sensitive to a change 
of the parameter a and therefore the corresponding lines cannot be clearly distinguished 
in the spectrum. 

For the consideration of the second family of lines we introduce new coordinates u = 
(x — y)/2 and v = [x -\- j/)/2, and consider a region that is bounded by the lines u = 0, 
M = I, u = and V = 1 and two circular arcs as shown in figure 7. Let us denote this region 
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by D. It contains four copies of the fundamental region of the biUiard which is indicated 
by thick hues. The second family of bouncing ball orbits runs along lines of constant u or 
constant v. In order that a wave function in the domain D is also an eigenfunction of the 
two considered symmetry classes of the fundamental region, it has to satisfy Dirichlet or 
Neumann boundary conditions, respectively, on the straight sections of the boundary of D 
and additionally it has to vanish on the diagonals v = u and v = 1 — u. 

We consider now a wave function which is mainly concentrated along the family of 
bouncing ball orbits that is parallel to the u-axis. This region is indicated by the tilted 
rectangular region in figure 7. Again we try to approximate the energy of this scarred 
wave function by an eigenvalue of the rectangular billiard. Depending on the considered 
symmetry class, the eigenfunction /(u, v) of the rectangular billiard has to satisfy Dirichlet 
or Neumann boundary conditions on its short sides and Dirichlet boundary conditions 
on its long sides (we consider the case a = 0). From f[u^v) one can obtain a wave 
function tl^{u^ v) in the fundamental region, i. e. one which vanishes on the lines v = u and 
u = 1 — M, by defining tl^{u^ v) = f{u^ v) — /(u, u) and requiring that f[l/2-\-z^l/2 — z) = 
/(1/2 — 2;, 1/2 + z). A consequence of the last condition is that the function /(u, v) must 
have the same parity with respect to refiection of its first coordinate at u = 1/2 as with 
respect to refiection of the second coordinate at v = 1/2. 




Figure 7: The (x, y)- and (u, u)-coordinate systems for a description of the two families of 
bouncing ball orbits in the quarter Sinai billiard with R = 0.5 and L = 2. The fundamental 
region of the Sinai billiard is drawn with thick lines. 

We now consider the odd symmetry class of the quarter Sinai billiard. For this case the 
function /(u, v) is given by /(u, v) = Cn^m sin[7rnu] sin[7rm(M — r)/(l — 2r)], where Cn^m is a 
normalization constant and r = i?/\/2 is the radius of the circle in the (u, u)-coordinates. 
The corresponding wave numbers are kn^m = VTA/n^/a^ + m^/6^/\/2 where a = 1 and 
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6 = (1 — 2r) and the factor -y^ has its origin in the Jacobian determinant of the coordinate 
transformation. The symmetry condition that /(u, v) has to satisfy restricts the possible 
values of m and n. It is satisfied if either m and n are both even or if they are both odd. 
As a consequence, the lowest excitation m = 1 in u-direction is only possible, when n is an 
odd number. In the left graph of figure 6 the wave numbers kn^m are indicated by empty 
triangles for the values n = 87 and 89 and m = 1 and one can see close to these positions 
horizontal lines in the spectrum. The positions of the lines are actually slightly below the 
values of kn^rm possibly because the wave functions spread slightly beyond the rectangular 
region. The lines are less pronounced than the lines of the first family and already for the 
cases m = 2, which correspond to even values of n, there are no clear lines in the spectrum. 

For the even symmetry class of the quarter Sinai billiard, the function /(u, v) is given 
by /(u, v) = Cn^m cos[7rnu] sin[7rm(M — r)/(l — 2r)]. Now the symmetry condition for /(u, v) 
requires that m is even if n is odd and vice versa. In the right graph of figure 6 the lines 
of the second family for m = 1 have the same qualitative properties as in the left graph, 
with the difference that they now appear for even values of n. The empty triangles in the 
right graph mark the values of kn^m for the values n = 86 and 88 and m = 1. 

Up to now, we considered only the position of the horizontal lines in the spectrum near 
the parameter value a = 0. We now discuss the fact that the lines are not strictly horizontal 
but that they drop slightly when a is increased. In order to describe the lines more 
accurately over the whole a-range we again use the simple model of the two rectangular 
billiards in figure 7. This time we impose mixed boundary conditions on those sides of 
the two rectangles that are tangent to the circular arcs in figure 7. The resulting wave 
numbers kn^m{ct) (for the appropriate values of n and m) have the following properties. 
They are monotonously decreasing functions of a, yet the decrease is very slight except 
near a = 7r/2, where it becomes stronger. At a = 7r/2 the lines have a distinctive negative 
curvature and negative slope. The model of the rectangular billiard is not good enough 
to give an exact description of the actual position of the lines in the spectrum, but the 
qualitative behavior of the functions kn^m{ct) that is described above can also be seen in 
the spectrum. 

5.2 The Velocity and Curvature Distributions 

We now come back to the distribution of velocities and curvatures and consider first the 
results for the parameter value a = 0. The velocity distribution is shown in figure 8 and 
is compared to a Gaussian of the same mean value and standard deviation. One notices 
that the distribution is not well approximated by the Gaussian curve. It is asymmetric and 
it has an excess of large velocities. This deviation from the random matrix result can be 
explained by the presence of the horizontal lines in figure 6 (see [12]). When the spectrum 
is unfolded, these lines go over into lines that have a large slope (except near a = 7r/2 as 
will be discussed later). They thus contribute to large positive velocities. Similar results 
were obtained for the hydrogen atom in a strong magnetic field, a system which has a large 
number of scarred states [12]. 

The curvature distribution that is shown in figure 9 is likewise not in agreement with 
the random matrix result. There is a large excess of small curvatures in comparison with 
the distribution of eq. (92). These results are again in qualitative agreements with results 
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Figure 8: The velocity distribution P{v) at a = in comparison with a Gaussian. 
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Figure 9: The curvature distribution P{c) at a = in comparison with the GOF- 
distribution of eq. (92). 
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for the hydrogen atom in a magnetic field and also with results for the stadium billiard 
where the parameter is the length of the straight section of the boundary [31, 12]. 

The high peak of the distribution at small curvatures can be explained by two effects. 
The horizontal lines in figure 6 contain long pieces where the curvature is almost zero. This 
is still the case after unfolding the levels and these parts of the spectrum thus contribute 
to an accumulation of small curvatures. Secondly, as has been noted above they are also 
responsible for an excess of large velocities. This excess leads to an increase of the variance 
of the velocity distribution and, by definition (91), to a decrease of the curvatures. 

Because of these effects it has been suggested that the excess of small curvatures can be 
taken as measure for the degree of scarring in a system [31, 12]. We will discuss this point 
in more detail when considering results for other values of the parameter a. 

It is possible that the small avoided level crossings along the horizontal lines lead also 
to a non-generic curvature distribution at large curvatures [12], but the number of levels 
has not been large enough in order to check this in detail. 

In order to see whether the discrepancy between the calculated distributions and the 
random matrix result is only due to levels in the vicinity of the horizontal lines, we disre- 
garded these levels and determined the distributions with the remaining levels. The levels 
which are taken out, however, have on average large velocities and therefore the average ve- 
locity of the remaining levels is different from zero. This effect can be removed by choosing 
a different unfolding mechanism which takes into account the fiuctuations in the spectrum 
due to the bouncing ball orbits. This is done by defining 



where Nbb[k) is the semiclassical contribution of the bouncing ball orbits to the spectral 
staircase. We would like to note that this unfolding procedure alone already leads to a more 
symmetric velocity distribution and reduces slightly the peak of the curvature distribution 
at small curvatures. 

We varied the number of levels that were removed and found that we had to neglect 
about 25 % of the spectrum until most of the non-generic contributions to the distributions 
were removed. The resulting velocity and curvature distributions, obtained by unfolding 
with (93), are shown in figures 10 and 11. The agreement with the random matrix results 
is now much better than before. One can still see some small deviations, but the main 
part of the non-generic contributions are removed. This indicates that wave functions that 
are scarred along the bouncing ball orbits are really the reason for the deviation from the 
random matrix results. 

We have calculated the velocity and curvature distributions also for several other values 
of the parameter a. Over a large range of this parameter the results are very similar to 
those for a = 0. The velocity distributions are asymmetric and have a tail at large velocities 
which is larger than for a Gaussian and the curvature distributions show an excess of small 
curvatures. Also, if one removes again the part of the spectrum which is mostly affected 
by the bouncing ball orbits and performs the unfolding with eq. (93) the agreement with 
the random matrix results becomes quite good. However, as one approaches the value 
a = 7r/2 the distributions become qualitatively different. In figure 12 the distribution of 
all velocities in the range 100 < < 300, unfolded without bouncing ball terms, is shown 
for a = 7r/2. It is compared to a Gaussian with the same mean and standard deviation. 
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Figure 10: The velocity distribution P{v) of a part ol the spectrum at a = in comparison 



with a Gaussian, as described in the text. 




Figure 11: The curvature distribution P{c) of a part of the spectrum at a = in comparison 
with the GOF-distribution of eq. (92), as described in the text. 
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Here the agreement with the Gaussian is better than for the distribution for a = in 
figure 8. The distribution is more symmetric, the peak at u = is lower and there is a 
smaller excess of large velocities. The curvature distribution is shown in figure 13. It differs 
strongly from the a = case. There are much less contributions from small curvatures 
and the deviations from the random matrix curve are mainly due to an asymmetry in the 
distribution. 

The different behavior near a = 7r/2 can again be explained by the apparent lines in 
the spectrum. As was discussed in the last section, these lines are not quite horizontal 
near a = 7r/2, but they have a negative slope and negative curvature. As a consequence, 
when the levels are unfolded they have smaller velocities and larger (negative) curvatures 
than in the case a = 0. This is exactly the difference that is observed in the velocity and 
curvature distributions. 

We also tried for a = 7r/2 to improve the agreement with the random matrix results 
by removing a part of the spectrum and unfolding with eq. (93). However, in this case the 
agreement does not improve. The peak of the velocity distribution goes down, but the 
distribution is slightly asymmetric and the peak of the curvature distribution decreases 
below that of the random matrix result. Presently, we do not have a satisfactory expla- 
nation for this result. Possible explanations are that due to the stronger a-dependence of 
the horizontal lines they have an infiuence also on more distant energy levels, or that the 
mixed boundary conditions do not behave generically near Neumann boundary conditions, 
at least in the considered energy range. 

A definite interpretation of our results would also require an examination of wave func- 
tions. However, in case that there are many scarred wave functions at a = 7r/2, our results 
seem to indicate, that scarred states do not necessarily have to lead to an accumulation 
of small curvatures. A prerequisite for such an accumulation is an almost linear change of 
the corresponding energy levels with the parameter a, which is not necessarily the case. 

Finally, we add a remark on the choice of the parameter a. The curvature distribution 
is not independent on this choice if it is evaluated within a restricted energy range. We 
chose a parameter dependence that is symmetric with respect to Dirichlet and Neumann 
boundary conditions, but in principle, one could define a different parameter A = f{a) 
and then the curvature distribution would change. For example we calculated at a = 7r/4 
the curvature distribution also for the K-dependent boundary conditions of eq. (4) and 
we got a slightly different result than for the a-dependence. (At a = 7r/2 there is no 
difference, since the Taylor expansion of the sine has no quadratic term.) Fxaminations of 
the universality of the curvature distribution thus require also a careful definition of the 
parameter dependence. 

6 Discussion 

In this paper we have extended the semiclassical methods for treating billiard systems to 
situations where the boundary conditions interpolate between the standard Dirichlet and 
Neumann conditions. The parametric dependence on the boundary condition was shown 
to be a very useful tool in the spectral analysis and in the context of the parametric 
spectral statistics, the variation of the boundary condition is similar to the variation of 
any external parameter in a Hamiltonian system. It should be born in mind, however. 
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V 

Figure 12: The velocity distribution P{v) at a = 7r/2 in comparison with a Gaussian. 




Figure 13: The curvature distribution P{c) at a 
distribution of eq. (92). 



7r/2 in comparison with the GOF- 
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that the boundary condition is strictly a wave (quantum) feature and as such, has no 
classical analogue. In this sense the system we consider here is different from the usual 
parameter dependent Hamiltonian systems, where the parameter has similar meanings in 
the quantum and classical representations. We shall address this point in the present 
section and show that, in one limit, there exists an alternative interpretation of the mixed 
boundary conditions. We will discuss also some results which can further illuminate some 
features of the semiclassical approximation for mixed boundary conditions. 

We shall show first, that the variation of the spectral density in the vicinity of the 
Dirichlet boundary condition can be related to a variation of the boundary itself. The 
later certainly has a similar meaning in both the classical and the quantum pictures. We 
rewrite the mixed boundary conditions, eq. (4), as 



V'(r ) + ednipi f ) = , r G S 



where e 



respect to e near e 



e > 0. This defines a Green function, G{ 



r' r "; e) 



(94) 

and its derivative with 



0"*" (Dirichlet case), g{r\r") 



can be easily shown to satisfy 



g{f ,f )+ ds{p) Go{f ,p) dng{p,f' 



ds{p) a^G'o(r ', p) dnGip, f"; e = 0) , (95) 



where G'o(r, r*') is the free Green function. Let us now define a uniform infiation of the 
boundary, which is affected by shifting each point p on the boundary a distance p along 
the normal: 

E(//) : p^ p + pn{p) , S . (96) 

Note, that this infiation is not just a rescaling of the coordinates , except for special cases 
like the circle billiard and regular polygons. Hence, the structure of the periodic orbits will 
in general change. Keeping the boundary condition to be Dirichlet (e = 0) for all values of 
one can show after some calculation, that dG{r\ r*"; v) j dv\i,=o satisfies the same integral 
equation (95) as for g{r\r"). Thus we establish the relation between Green functions 



dG{f', f"- e, v) 



de 



dG{f', f"- e, v) 



t=u=0 



dv 



t=u=0 



To derive the relation between the spectral densities we use for A; > 

2k 



Zk f 

d{k\t,v) = lim / dr Im Gfr, r*; A;^ + 

^ ^ TT 5^0+ J A 



(97) 



(98) 



The derivative with respect to e affects only the integrand, while the derivative with respect 
to V affects also the integration domain. However, since we consider Dirichlet boundary 
conditions (e = 0), we have G{f^f; e = 0) = for r*G S, and the second contribution is 0. 
Therefore we conclude: 



dd{k; e, u) 



de 



dd{k; e, u) 



e=iy=0 



dv 



(99) 



e=iy=0 



This exact result becomes intuitively clear in the semiclassical limit: For the smooth part 
of the level density we can use dA/ dp\,y=o = Z, where A, Z are the billiards area and 
perimeter, respectively. When this is inserted into eq. (7) we get to leading order: 



dd{k; e, u) 
dl 



e=iy=0 



dd{k; e, u) 
dv 



e=iy=0 



Lk 
27 



(100) 
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(For smooth billiards one gets also the next term by similar methods.) As for dosc^ con- 
sider only the contributions of generic periodic orbits which are far from tangency to the 
boundary. In the Gutzwiller sum, only the fast varying derivatives of the actions kL^ 
with respect to p will be important in the semiclassical limit. Simple geometry shows that 
dL^ I dv\u=o = 2J2iZi cos^!^, and together with the relation d/de\^=o = b d / da\a=o one gets 
the desired equivalence (see eq. (86) ). 

Even though the boundary condition is a pure quantum feature, its expression in the 
semiclassical spectral density depends on the underlying classical dynamics and deserves 
some more discussion. In particular, we shall study now the phase factor which multiplies 
the contributions of the generic unstable periodic orbits (9). The main feature of this phase 
is that for the mixed boundary conditions it depends on both, k and the angles of incidence 
Oi^i = 1, . . . n along the periodic orbit. This should be contrasted with the Dirichlet and 
Neumann boundary conditions, where a reflection leads to a simple factor of ( — 1) or (+1), 
respectively. 

Let us consider the phase factor for very long periodic orbits. We are interested in orbits 
that fill the phase space uniformly, so that the phase factor can be calculated approximately 
by an ergodic averaging over phase space. The appropriate coordinates for doing this are 
the length along the perimeter s and the canonical momentum p = cos where (p is the 
angle of the refiected trajectory with the (directed) tangent. The phase $ is then obtained 
as 

fk \ 

$ = 2 ^ arctan ( — cos 9i 1 

1 / 
dp arctan ■ 
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The number of refiections on the boundary n can be expressed byn = // <c>, where 
/ is the length of the orbit and < c > is the mean chord length, i. e. the mean distance 
between two consecutive refiections on the boundary. Denoting by c(s,(^) the length of a 
chord that starts at s with an angle (p with respect to the tangent, the mean chord length 
is obtained by 

<c> = — - [ ds [ dp s'm Lp c{s,Lp) 
2L Jo Jo 

I I-2TT l-L 

2L 
_ ttA 

~L 

where a transformation to an angle a with respect to a fixed direction in space has been 
made and c(s, p) = for directions that point outside of the billiard. The integral over ds 
gave the area A, since ds • sin p is the component of ds orthogonal to the chord. With this 
result the phase $ is given by 



/ da / ds sin (^(a) c(s, (^(a) ) 
Jo Jo 

(102) 
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This additional phase can be interpreted as a small correction to the classical action / • k 
which is the leading term in the semiclassical phase. The expression which multiplies / in 
eq. (103) leads to a mean shift in the spectrum with respect to the Dirichlet value: 

On the other hand the mean spectral shift due to a change in the mean number of levels 
is semiclassically given to leading order by (7) 

kn - C = -[Nikn) - N^ikn)]Ak = -j (^^1 + (^) ^ - , (105) 

where Ak = l/d[k) is the mean distance between wave numbers. The result is identical 
to eq. (104). It shows that the ergodic average over the additional phase is responsible 
for the mean shift of the energy levels with respect to the Dirichlet case. The deviations 
of the phase factors of periodic orbits from the ergodic mean thus lead to fluctuations of 
the energy levels around their mean shift. The fact that the mean shift in the energy 
levels, which is derived from considering very short orbits, can also be obtained from the 
ergodic average of very long orbits is another example of the bootstrapping property in 
periodic-orbit theory [36]. 

The last comment is a word of caution. The mixed boundary conditions were shown 
to appear in the expressions for the smooth and the oscillatory components of the spectral 
density. The theory presented here for the smooth part gives the first few leading terms 
in the asymptotic expansion of d[k) in powers of k~^ . This is not true for the oscillatory 
contribution, where only the leading contribution was considered. 
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A The Smooth Part of the Spectral Density for the 
Semi-circle Billiard 

In this appendix, the smooth part of the spectral staircase N[k) is derived for a semi- 
circle with mixed boundary conditions on the circular part of the boundary and Dirichlet 
or Neumann boundary conditions on the diameter, respectively. The derivation is done 
by using the methods of section 2.1. Erom the results the contribution of a 90° corner is 
extracted. 
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The Green function for the semi-circle billiard is obtained from that of the circle billiard 
by the method of images: 

G{f, r s^) = G%f, r s^) T G%f, Rr s^) , (106) 

where the upper and lower sign correspond to Dirichlet and Neumann boundary conditions 
on the diameter, respectively and the operator R denotes the reflection on the line y = 0. 
The upper index c will from now on denote the quantities for the full circle billiard. From 
eq. (106) the Green function for the semi-circle follows as 

-1 oo 

G{f,f',s') = -— J2 In{sr<)[K.,{sry) + aJ.,{sry)][cosn{e-e')Tcosn{e + e')], (107) 

n=-oo 

Subtracting the free Green function and taking the trace results in 

Jo Jo 

= ^K^{s') T ^[/(O, s) - Io{sR)Ko{sR) + I^{sR)K^{sR)] . (108) 

Now, the asymptotic expressions for the Bessel functions for large arguments are in- 
serted and after expanding in powers of l/[sR) one obtains 



2 



k/s-1 1 



2{sRy {k/s + 1) sR 



(109) 



From this the contributions to the mean level density and to the mean spectral staircase 
follow as 

d(k) = -d%k)T — T^ ^ — (110) 

^ ' 2 ^ ^ ^ 27r ^ 27rA; ^ ^ (^^^y ^ ' 

and 

1 - A- /? 1 A- 1 
iV(A;) = -iV=(A;) T — T (— arctan ), (111) 

where the constant term has been determined as discussed in section 2.1. The term 
=pA;i?/(27r) in eq. (Ill) is the contribution of the diameter to N[k). The last term in 
eq. (Ill) is the contribution of the two corners and thus it follows that the contribution 
of a 90° corner with mixed boundary conditions on one side and Neumann and Dirichlet 
boundary conditions, respectively, on the other side is given by 

=F(— arctan— —) . (112) 



References 

[1] N. S. Koshlyakov, M. M. Smirnov, and F. B. Gliner. Differential Equations of Math- 
ematical Physics. North-Holland, Amsterdam, 1964. 

[2] P. M. Morse and H. Feshbach. Methods of Theoretical Physics, Part II. McGraw - 
Hill, New York, 1953. Chapter 9.2. 



43 



[3] R. Balian and C. Bloch. Distribution of eigenfrequencies for the wave equation in a 
finite domain: I. Three-dimensional problem with smooth boundary surface. Ann. 
Phys., 60:40f-447, f970. Errata in Ann. Phys., 84:559-563, f974. 

[4] M. C. Gutzwiller. Periodic orbits and classical quantization conditions. J. Math. 
Phys., f2:343-358, f97f. 

[5] M. C. Gutzwiller. Chaos in Classical and Quantum Mechanics. Springer, New York, 
f990. 

[6] M. V. Berry and M. Tabor. Closed orbits and the regular bound spectrum. Proc. R. 
Soc. Lond. A., 349:f0f-f23, f976. 

[7] M. V. Berry and M. Tabor. Calculating the bound spectrum by path summation in 
action-angle variables. J. Phys. A: Math. Gen., f0:37f-379, f977. 

[8] M. V. Berry. Quantizing a classically ergodic system: Sinai billiard and the KKR 
method. Ann. Phys., f3f:f63-2f6, f98f. 

[9] M. Sieber, U. Smilansky, S. C. Creagh, and R. G. Littlejohn. Non-generic spectral 
statistics in the quantized stadium billiard. J. Phys. A: Math. Gen., 26:62f 7-6230, 
f993. 

[fO] P. Gaspard, S. A. Rice, H. J. Mikeska, and K. Nakamura. Parametric motion of energy 
levels: Curvature distribution. Phys. Rev. A, 42:40f 5-4027, f990. 

[ff] J. Goldberg, U. Smilansky, M. V. Berry, W. Schweizer, G. Wunner, and G. Zeller. 
The parametric number variance. Nonlinearity, 4:f-f4, f99f . 

[f 2] J. Zakrzewski and D. Delande. Parametric motion of energy levels in quantum chaotic 
systems. I. Curvature distributions. Phys. Rev. E, 47:f 650-f 664, f993. 

[f3] B. D. Simons and B. L. Altshuler. Universalities in the spectra of disordered and 
chaotic systems. Phys. Rev. B, 48:5422-5438, f993. 

[f4] K. Stewartson and R. T. Waechter. On hearing the shape of a drum: Further results. 
Proc. Camb. Phil. Soc, 69:353-363, f97f. 

[f5] B. Dietz and U. Smilansky. A scattering approach to the quantization of billiards - 
the inside-outside duality. Chaos, 3:58f-590, f993. 

[f6] B. Dietz and U. Smilansky. Scattering from a square obstacle, preprint, f993. 

[f7] J. -P. Eckmann and C.-A. Pillet. Spectral duality for planar billiards. Commun. Math. 
Phys., 1995. in press. 

[f8] W. Magnus, E. Oberhettinger, and R. P. Soni. Formulas and Theorems for the Special 
Functions of Mathematical Physics. Springer, Berlin, f966. 

[f 9] M. V. Berry and C. J. Howls. High orders of the Weyl expansion for quantum billiards: 
Resurgence of periodic orbits, and the Stokes phenomenon, preprint, f994. 



44 



[20] M. Abramowitz and I. A. Segun, editors. Handbook of Mathematical Functions. Dover, 
New York, 1965. 

[21] H. P. Baltes and E. R. Hilf. Spectra of Finite Systems. Bibliographisches Institut, 
Mannheim, 1976. 

[22] E. Doron and U. Smilansky. Semiclassical quantization of chaotic biUiards: a scattering 
theory approach. Nonlinearity, 5:1055-1084, 1992. 

[23] U. Smilansky and I. Ussishkin. in preparation. 

[24] J. Korringa. On the calculation of the energy of a Bloch wave in a metal. Physica, 
13:393-400, 1947. 

[25] W. Kohn and N. Rostoker. Solution of the Schrodinger equation in periodic lattices 
with an application to metallic Lithium. Phys. Rev., 94:1111-1120, 1954. 

[26] E. S. Ham and B. Segall. Energy bands in periodic lattices - Green's function method. 
Physical Review, 124(6):1786-1796, 1961. 

[27] H. Primack and U. Smilansky. Quantal consequences of perturbations which destroy 
structurally unstable orbits in chaotic billiards. J. Phys. A, 27:4439-4448, 1994. 

[28] H. Primack and U. Smilansky. Quantization of the three dimensional Sinai billiard, 
submitted for publication, 1995. 

[29] H. Schanz and U. Smilansky. Quantization of Sinai's billiard - A scattering approach, 
in press, 1995. 

[30] D. Saher, E. Haake, and P. Gaspard. Parametric-curvature distribution in quantum 
kicked tops. Phys. Rev. A, 44:7841-7843, 1991. 

[31] T. Takami and H. Hasegawa. Curvature distribution of chaotic quantum systems: 
Universality and nonuniversality. Phys. Rev. Lett., 68:419-422, 1992. 

[32] E. von Oppen. Exact distribution of eigenvalue curvatures of chaotic quantum systems. 
Phys. Rev. Lett, 73:798-801, 1994. 

[33] Y. V. Eyodorov and H.-J. Sommers. Universality of "level curvature" distribution for 
large random matrices: Systematic analytical approaches, preprint, 1994. 

[34] K. Nakamura. Quantum Chaos. Cambridge University Press, 1993. 

[35] M. KoUmann, J. Stein, U. Stoffregen, H.-J. Stockmann, and B. Eckhardt. Periodic 
orbit analysis of billiard level dynamics. Phys. Rev. E, 49:R1-R4, 1994. 

[36] M. Berry. Some quantum-to-classical asymptotics. In M. J. Giannoni, A. Voros, and 
J. Zinn- Justin, editors, Les Houches 1989 Session LII, Chaos and Quantum Physics, 
pages 250-303. North-Holland, Amsterdam, 1991. 



45 



